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APPENDIX A 

The following is an alternative implementation of the present invention. The 
following code is neither annotated nor further described, but is readily understood by 
ordinarily skilled software developers, especially given the thorough description the invention 
in previous subsections. The following code contains material which is subject to copyright 
protection. The copyright owner has no objection to the facsimile reproduction by anyone of 
the patent document or the patent disclosure, as it appears in the Patent and Trademark Office 
patent file or records, but otherwise reserves all copyright rights whatsoever. 



#define UPDATE_CENTROID 
#define USE_MORPH_CLOSE 

#define USE_SCALAR_KMEANS_INITIALIZATION 
#define USE_LPF_IN_KMEANS 

# define MY_NEG_INF -HUGE_VAL/ 100 

#define MY_POS_INF HUGE_VAL/100 

#define ALAW_COMPAND_PARAM 25 // parameter A for Alaw 

Companding . 

#define MAX_INTENSITy_VALUE 65535 // input signal peak magnitude 

inline double 

log_gaussian_prob (double x, double mean, double istd) 

{ //un-normalized and unsealed by std, dev. log of gaussian prob 

if (mean==0) return MY_NEG_INF; 

else if (istd==0) return 0.0; 

else { 

double z= (x-mean) *istd; 
return {-z*z) ; 

} 

> 

inline double 

pos_prob(int type, double r, double nr) 

{ //un-normalized position prob for FG (type=l) or BG (type=0) 
const double RFACTOR=0.5; 

static double nr_L0_N0MR, nr_HI_NOMR, nr_DI_NOMR; 

switch (type) { 
case 0 : 

i f ( r<=nr_LO_NOMR ) return 0.0; 

else if {r<=nr_HI_NOMR) return (r-nr_LO_NOMR) * {nr_DI_NOMR) ; 
else return 1.0; 
case 1: 

if (r<=nr_LO_NOMR) return 1.0; 

else if (r<=nr_HI_NOMR) return (nr_HI_NOMR-r) * (nr_DI_NOMR) ; 
else return 0.0; 
case 2 : 

nr_LO_NOMR=nr* ( 1 . 0-RFACTOR) ; 
nr_HI_NOMR=nr* ( 1 . 0+RFACTOR) ; 
nr_DI_NOMR=l .0 / (nr_HI_NOMR-nr_LO_NOMR) ; 
return 0.0; 

> 

return 0.0; 
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} 

inline double 

pos_probCu ( int type, double r, double nr) 
5 { //un-normalized position prob for FG (type=l) or BG (type=0) 
const double RFACTOR=0 . 5 ; 

static double nr_LO_NOMR, nr_HI_NOMR, inult_f actor ; 

switch (type) { 
10 case 0: 

if (r<=nr__LO_NOMR) return 0.0; 
else if (r<=nr_HI_NOMR) { 
double t=(r-nr); 
return 0 . 5-inult_f actor*t*t*t ; 

15 } 

else return 1.0; 
case 1 : 

if (r<=nr_LO_NOMR) return 1.0; 
else if (r<=nr_HI_NOMR) { 
20 double t=(r-nr); 

return 0. 5+inult_f actor*t*t*t ; 

} 

else return 0.0; 
case 2 : 

25 mult^f actor= (nr*RFACTOR) ; 

mult_f actor=0 . 5/ (mult_factor*mult_factor*mult_f actor) ; 
nr_LO_NOMR=nr* ( 1 . 0-RFACTOR) ; 
nr_HI_NOMR=nr* ( 1 . 0+RFACTOR) ; 
return 0.0; 

30 } 

return 0.0; 

} 

void 

35 GetlnitialClassif icationNomRaddong iXStart, long iXStop, long lYStart, 
long lYStop, 

double dCentroid_x, double dCentroid_y, double nom_rad, 
CTArray<bool> &pBin) 
{ 

40 long p,x,y; 

double r; 

for (p=0 , y=lYStart ; y<=lYS top ; y++ ) 

for (x=lXStart;x<=lXStop;x++,p++) { 

r=sqrt ( (x-dCentroid__x) * (x-dCentroid_x) + (y- 
45 dCentroid_:y) * (y-dCentroid_y ) ) ; 

pBin[p] = (r<=nom_rad) ; 

} 

} 

50 void 

DoHoleFill {long iXStart, long IXStop, long lYStart, long lYStop, 

CTArray<bool> &pBin) 

{ 

long lRowSize=lXStop-lXStart+l; 
55 long lNuinPixels=pBin . GetSize ( ) ; 

long p,x,y; 
#ifdef USE_MORPH_CLOSE 

// Morphological Close operation in 3x3 cross window 

CTArray<bool> tBin; 
60 tBin. Set Size (INumPixels) ; 

// Copy edge pixels 

for (p=0;p<lRowSize;p++) tBin [p] =pBin [p] ; 
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for (y=lYStart +1 ; y<lYStop ; y++ , p+=lRowSize ) { 

tBin[p] =pBin[p] ; tBin [p+lRowSize-1 ] =pBin [p+lRowSize-l ] ; 

} 

for ( ;p<lN\amPixels;p++) tBin[p] =pBin [p] ; 
5 // Dilate 

for (p=lRowSize+l, y=lYStart+l;y<lYStop;y++,p+=2 ) 
for (x=lXStart+l ; x<lXStop ; x++ , p++ ) 

tBin[p] = (pBin[p] | | pBin[p-l] | | pBin[p+l] | | pBin[p- 
IRowSize] I I pBin[p+lRowSizel ) ; 
10 // Erode 

for (p=lRowSi ze+1 , y=lYStart+l ; y<lYStop ; y++ , p+=2 ) 
for (x=lXStart+l ; x<lXStop; x++ , p++ ) 

pBin[p] = (tBin[p] && tBin[p-l] && tBin[p+l) && tBin[p- 
IRowSize] && tBin[p+lRowSize] ) ; 
15 #else 

// A simple filling operation to fill stray BG pixels surrounded by 
// 4-connected FG pixels and vice-versa 
for (p=lRowSi ze+1 , y=lYStart + l ; y<lYStop ; y++ , p+=2 ) { 
for (x=lXStart+l;x<lXStop;x++,p++) { 
20 if {pBin[p]==false) 

pBin[p] = {pBin[p-l] && pBin[p+l] && pBin[p- 
IRowSize] && pBin[p+lRowSize] ) ; 

else 

pBin[p] = {pBin[p-l] | | pBin[p+l] | | pBin[p- 
25 IRowSize] II pBin[p+lRowSize] ) ; 
} 

> 

#endif 

return; 

30 } 

// param provides the rau or A value for companding. InPk provides the input 

signal peak magnitude. 

bool 

35 ALawCompand ( CTArray< CTArray<unsigned short> > &ppPixelData, int param, 
unsigned short InPk) 
{ 

double InApl, InPkdA; 
long p, c; 

40 // CTArray< CTArray<unsigned short> > ppPixelDataALaw; 
long lNumColors=ppPixelData . GetSize ( ) ; 
long lNumPixels=ppPixelData [0 ] .GetSizeO; 
InApl = log (param) + 1; 
InPkdA = InPk / param; 
45 // A-law compressor 

for (p=0;p<lNumPixels;p++) 

for (c=0;c<lNumColors;c++) 
{ 

if ( (ppPixelData[c] [p] ) <= InPkdA ) 
50 ppPixelData [c] [p]= (unsigned short) (param / InApl * 

ppPixelData[c] [p]+0.5); 

else if ( (ppPixelData [c] [p] ) > InPkdA) ; 

ppPixelData [c] [p]= (unsigned short) (InPk / InApl * 
(1 + ppPixelData [c] [p] ) / InPkdA + 0.5); 
55 } 

return 1 ; 

} 

bool 

GetlnitialClassif icationScalarKMeans (long iXStart, long IXStop, long 
60 lYStart, long lYStop, 

CTArray< 

CTArray<unsigned short> > &ppPixelData, CTArray<bool> &pBin, 
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int 

nuin_no i s e_ t hr e s h ) 
{ 

const long MAX_KM_ITER=32; 
5 const long MAX_NO_TRY=3 ; 

const double CONVERGENCE_FACTOR=0 . 001; 
const long weakthresh=20 ; 

long lNuiciColors=ppPixelData .GetSize ( ) ; 
10 long lNiainPixels=ppPixelData [0] .GetSize () ; 

unsigned long 1[2], lsuin[2] , It, Itprev; 

long p, c,n[2] ,m[2] , iNumlter, INumTry; 

bool bIsConverged, bIsDone , ex; 

CTArray<unsigned long> pPixelSumData; 
15 CTArray<bool> mBin; 

pPixelSioinData . SetSize ( INumPixels ) ; 
mBin. Set Size (INumPixels) ; 

20 //Computes sum of pixels for 2 -means 

for (p=0;p<lNumPixels;p++) 

for (pPixelSumData [p] =0 , c=0 ; c<lNumColors ; C++ ) 

pPixelSumData [p] += ( (unsigned long) ppPixelData [c] [p] ) ; 

25 #ifdef USE_IiPF_IN_KMEANS 

//Low pass filter data 

CTArray<unsigned long> pPixelS\xmFiltData; 
pPixelSumFiltData . SetSize { INumPixels) ; 
long lRowSize=lXStop-lXStart+l , x, y ; 
30 for {p=0 ;p<lRowSize;p++) pPixelSumFiltData [p] =pPixelS\ainData [p] ; 

for (y=lYStart+l;y<lYStop;y++,p+=lRowSize) { 
pPixelSumFiltData [p] =pPixelSumData [p] ; 
pPixelSumFiltData [p+lRowSize-1] =pPixelS\imData [p+lRowSize-1 ] ; 
} 

35 for ( ;p< INumPixels ;p++) pPixelSumFiltData [p] =pPixelSumData [p] ; 

for (p=lRowSize+l,y=lYStart+l;y<lYStop;y++,p+=2) { 
for {x=lXStart+l;x<lXStop;x++,p++) { 

pPixelSumFiltData [p] = ( { (pPixelSumData [p] «2 ) +pPixelSumData [p+1] +pPixe 
40 lSumData[p-l]+ 

pPixelSumData [p+lRowSize] +pPixelSumData [p-lRowSize] +4 ) >>3 ) ; 
} 

} 

45 for (p=0;p<lNuinPixels;p++) pPixelSumData [p] = (pPixelSumFiltData [p] ) ; 

#endif 

for (p=0 ;p<lNumPixels;p++) mBin[p]=l; 

50 bIsDone=f alse; lNximTry=0; 

while (!bIsDone ) { 

//Computes max and min of sum 
l[0] = (6553 5*lNumColors) ;1[1] = (0) ; 
for (p=0 ;p<lNumPixels ;p++ ) { 
55 if (mBin[p]) { 

1 [1] = (pPixelSumData [p3>l [1] ?pPixelSumData [p] :1 [1] ) ; 
1 EO] = (pPixelSumData [p]<l [0] ?pPixelSumData [p] :1 [0] ) ; 

} 

} 

60 if (ltO)==l[l]) return 0; 

//Do k-means for 2 classes 
lt=(l[0]+l[l]) ; 
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bIsConverged=f alse ; lNuinIter=0 ; 

while ( IblsConverged && lNuinIter<MAX_KM_ITER) { 
ltprev=lt ; 

lsum[0] =lsuin[l)=0 . 0; 

n[0]=n[l]=0; 

m[0]=m[l]=0; 

for (p=0 ;p<lNumPixels;p++) { 

pBin[p] = ( (pPixelSumData [p] «1) >lt) ; 
n(pBin[p] ]++; 
if (mBin[p]) { 

m[pBin[p] ] ++; 

lsuin[pBin[p] ] + = (double) pPixelSiunData [p] ; 

} 

} 

if (m[0]) 1[0]=( (double) lsum[0] /m[0]+0 .5) ; 
if (m[l]) 1[1] = ( (double) lsum[l)/in[ 11+0.5) ; 
lt=(l[0]+l[l] ) ; 

bIsConverged= ( f abs ( (double) It/ltprev- 
1) <CONVERGENCE_F ACTOR) ; 

lNumIter++; 

} 

if (bIsConverged==0 ) 

print f ( "DebugXn" ) ; 
if (n [ 1] <nuin_noise_thresh || n [ 0 ] <num_noise_thresli) { 

ex= (n [ 1 ] <nuin_noise_ thresh) ; 

for (p=0 ;p<lNumPixels;p++) 
inBin[p] = (mBin [p] ? (pBin[p] "^ex) : 0 ) ; 
} 

else 

bIsDone=true ; 
lNumTry++; 

} 

return (1 [1] -1 [01 > (weakthresh*lNuir\Colors ) ) ; 



long 

ComputeMask (CTArray< CTArray<unsigned short > > &ppPixelData, 

long IXStart, long IXStop, long lYStart, long lYStop, 

double noin_rad, 

double &dCentroid_x, double &dCentroid_y, bool 

&bIsSpotShifted) 
{ 

const long MAX_EM_ITER= 5 ; 

const long REJECT_ITER=3 ; 

const double REJECT_Z=3 . 0 ; 

const double REJECT_Z_STEP=2 . 0 ; 

const double MAX_REJECT_RATIO=0 . 75 ; 

const double KMEANS_RAD_RE JECT_FACTOR= 0 . 6 ; 

const double KMEANS_CEN_REJECT_FACTOR=l . 5 ; 

const double KMEANS_MOI_REJECT_FACTOR=2 . 0 ; 

const double ISSPOT_FACTOR=0. 75 ; 

const double ISWEAKSPOT_FACTOR=3 . 0 ; 

const double GRID_MOVED_REJECT_FACTOR=l . 0 ; 

static double zero=0; 

long p, X, y, c, lNumIter=0 , n [2 ] , nR [2 ] ; 
bool bIsConverged=f alse; 

long lNuinPixels= (lXStop-lXStart+1) * ( lYStop-lYStart+1 ) ; 
long lNuinColors=ppPixelData .Get Size ( ) ; 
double dRadius; 
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double gridCentroid_x, gridCentroid_y; 



gridCentroi d_x= dCen t r o i d_x ; 
gr i dCen t r o i d_y = dCent r o i d_y ; 



CTArray<bool> pBin, pBinPrev; 
CTArray<double> pPosProb; //Prob(Fg/r) 

CTArray< CTArray<double> > dAve, dStdDev, dSum, dSumSq, dSLimR, 
dSumSqR; 

pBin.SetSizedNumPixels) ; 
pBinPrev. Set Size (INumPixels) ; 

pPosProb. SetSize (iNumPixels) ; 

dAve . SetSize (2 ) ; 

dAve[0] . SetSize ( INumColors ) ; 

dAve[l] . SetSize (INumColors) ; 

dStdDev. SetSize (2) ; 

dStdDev[0] . SetSize ( INumColors ) ; 

dStdDev [1] . SetSize ( INumColors ) ; 

dSum. SetSize (2 ) ; 

dSum[0] . SetSize (INumColors) ; 

dSum[l] . SetSize (INumColors) ; 

dSumSq. SetSize (2 ) ; 

dSumSq[0] . SetSize ( INumColors ) ; 

dSumSq[l] . SetSize ( INumColors ) ; 

dSiunR . SetSize ( 2 ) ; 

dSumR[0] . SetSize (INumColors) ; 

dSumR[l] . SetSize (INumColors) ; 

dSumSqR. SetSize (2 ) ; 

dSumSqR[0] . SetSize (INumColors) ; 

dSiimSqR[l] . SetSize ( INumColors) ; 

ALawCompand(ppPixelData, ALAW_COMPAND_PARAM , 65535); 

#i f def USE_SCAIiAR_KMEANS_INITIALIZATION 

int minl= (int) (3 . 14159*pow( (1 . 0- 
KMEANS_RAD_REJECT_F ACTOR) *nom_rad, 2 . 0 ) +0 . 5 ) ; 

int minO = INumPixels - 
(int) (3 . 14159*pow( ( 1 . 0+KMEANS_RAD_REJECT_FACTOR/2 ) *nom_rad, 2 . 0 ) +0 . 5 ) ; 

minO= (minO<0?0 :minO) ; 

if ( IGetlnitialClassif icationScalarKMeans (IXStart , IXStop, lYStart, 
lYStop, ppPixelData, pBin, (minl<minO?minl iminO ) ) ) 
{ 

return 0 ; 

} 

DoHoleFill (IXStart , iXStop, lYStart, lYStop, pBin) ; 
for (n[0] =n[l] =0 , p=0 ;p<lNumPixels ;p++) n[pBin[p] ] ++; 
dRadius=sqrt ( (double ) n [ 1 ] /3 . 14159 ) ; 
bool acceptKM; 

if (dRadius> ( 1 . 0-KMRANS_RAD_REJECT_FACTOR) *nom_rad && 
dRadius< ( 1 . 0+KMEANS_RAD_REJECT_F ACTOR) *nom_rad) 
acceptKM=l; 

else 

{ 

return 0; 

} 
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double dtCentroid_x, dtCentroid y ; 
if (acceptKM) { 

//Compute centroid of fg region 

long suin_x=0, sum_y=0; 

double suinsq_r=0 . 0 ; 

for (p=0,y=lYStart;y<=lYStop;y++) { 

for (x=lXStart ;x<=lXStop;x++,p++) { 
if (pBin[p]) { 

s uni_x+ =x ; s uin_y + =y ; 

} 

} 

} 

if (n[l]) { 

dtCentroid_x = ( (double ) suin__x/n [ 1 ]) ; 

dtCentroid v = ( (double) sum__y/n [ 1] ) ; 

if (sqrt ( (dtCentroid_x-dCentroid__x) * (dtCentroid_x- 
dCentroid_x) + (dtCentroid y-dCentroid y) * (dtCentroid_y-dCentroid_y) ) > 
KMEANS_CEN_REJECT_FACTOR*nom_rad ) 

{ 

return 0 ; 
bIsSpotShif ted=l ; 

} 

sums<i_r=0 ; 

for {p=:0 , y=lYS tart ; y< = lYStop ; y++ ) { 

for (x=lXStart ;x<=lXStop;x++,p++) { 

if (pBin [p] ) suinsci_r+= {x-dtCentroid_x) * (x- 
dtCentroid_x) + (y-dtCentroid y) * (y-dtCentroid_y) ; 

} 

} 

if (sunisq„r*3 . 14159> KMEANS_MOI_REJECT_FACTOR*n 1 1 ] *n[l] ) 
{ 

return 0; 
bIsSpotShifted=l; 

} 

else acceptKM =1; 

> 

//Accept nom__rad and centroid based on K-Means 
noin_rad=dRadius ; 
dCentroid_x=dtCentroid_x ; 
dCentroid v=dtCentroid v; 

} 

else 

#endif 

GetlnitialClassif icationNomRaddXStart , IXStop, lYStart, lYStop, 

dCentroid_x, 

dCentroid nom_rad/ pBin) ; 

pos__prob (2,0, nom_rad) ; 

//Computes radii and normalization factors 
for (p=0 , y=lYStart ; y<=lYStop ; y++ ) 

for (x=lXStart ;x<=lXStop;x++ ,p++) { 

dRadius=sqrt ( (x-dCentroid_x) * {x-dCentroid_x) + (y- 
dCentroid v) * (y-dCentroid v) ) ; 

pPosProb [p] -pos__prob ( 1 , dRadius , nom_rad) ; 

} 



//EM iterations start 

while ( IblsConverged && lNumIter<MAX__EM_ITER ) { 
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double metric [2], precomp_ratio, reject_z(2]; 
double t ; 

//Computes averages and std. devs . of two classes 

5 

dSum[0] . Initialize (zero) ; 
dSumEl] . Initialize (zero) ; 
dSumSq(0] . Initialize ( zero) ; 
dS\iinSq[l] . Initialize ( zero) ; 
10 for (c=0;c<lNumColors; C++) { 

for (n[0] =n[l] =p=0 ;p<lNuinPixels;p++) { 
n[pBin[p] ] ++; 

dSuin[pBintp] ] [c ]+= (double) ppPixelData [c ] [p] ; 

15 dSiunSq[pBin[p] ] [c] += (double) ppPixelData [c] [p] * (double) ppPixelData [c] [ 

p]; 

} 

t = dAve[0] [c] = (n[0] ?dSuin[0] [c] /n[0] :0) ; 
t=dAve[l] [c] = (n[l]?dSuin[l] [c]/n[l] :0) ; 
20 t=dStdDev[0] [c]={n[Q]>l?sqrt (dSumSqlO] [c]- 

dAve[0] [c] *dAve[0] [c ] ) / (n [ 0 ] -1 ) : 0 ) ; 

t=dStdDev[l] [ c ] = (n [ 1] >l?sqrt (dSumSq[ 1] [c]- 
dAve[l] [c]*dAve[l] [c] ) / (n[l] -1) : 0 ) ; 
} 

25 

//reject outliers and recompute 
for (c=0 ; c<lNumColors ; C++) { 

reject_z [0] =REJECT_Z*dStdDev[0] [c] ; 
reject_z[l]=REJECT_Z*dStdDev[ll [c] ; ; 
30 for (int 1=0 ; 1<REJECT_ITER; 1++) { 

dSumREO] tc]=dSuin[0] [c] ; 
dSuinR[l] [c]=dSum[l] [c] ; 
dS\jmSqR[0] [c ] =dSuitiSq [ 0 ] [c] ; 
dSiimSqR[l} [ c ] =dS\ainSq [ 1 ] [c] ; 
35 for (nR[0}=nR[l] =p=0;p<lNuinPixels;p++) { 

if (fabs( ( (double) ppPixelData [c] [pl- 
dAve[pBin[p] ] [cl ) ) >=reject_z [pBin [p] ] ) { 

dSumR [pBin[p] ] [c] - 

= (double) ppPixelData [c] [p] ; 
40 dSumSqR[pBin[p] ] [c]- 

= (double) ppPixelData [c] [p] * (double) ppPixelData [c] [p] ; 

} else 

nR[pBin[p] ] ++; 

} 

45 if (nR[0]>n[0] *MAX_REJECT_RATIO) { 

t=dAve[0] Cc3=(nR[0] ?dSumRCO] [c]/nR[0] :dAve[0] [c] ) ; 

t=dStdDev[0] [c] = (nR[0] >l?sqrt (dS\HaSqR[0] [c]- 
dAve[0] [c] *dAvet0] [c] ) / (nR[0]-l) :dStdDev[0] [c3 ) ; 
50 } 

else reject_z[0]+=REJECT_Z_STEP*dStdDev[0] [c] ; 

if (nR[l]>n[l] *MAX_REJECT_RATIO) { 

55 t=dAve[l] [c]=(nR[l] ?dSumR[l] [c] /nR[l] :dAve[l] [c] ) ; 

t=dStdDev[l] [c]=(nR[13>l?sqrt(dSumSqR[l] [c]- 
dAve[l] [c] *dAveIl] [c] ) / (nR[ll-l) :dStdDev[l] [c] ) ; 

} 

else reject_z [l]+=REJECT_Z_STEP*dStdDev[l] [c] ; 

60 > 
} 
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bool bIsSpot=f alse; 

bool bIsStrongSpot=f alse; 

for (c = 0 ; c<lNuinColors ; C + +) { 

bIsSpot|=(fabs (dAve[l] [c]- 
dAve[0] [c] )>(dStdDev[l] [c] +dStdDev[0 ] [c] ) *ISSPOT_FACTOR) ; 

bIsStrongSpot I = (fabs (dAve [1] [c]- 
dAve[0] [C] )>(dStdDev[l] [c] +dStdDev[0] [c] ) *ISWEAKSPOT_FACTOR) ; 
} 

i f ( ! bIsSpot ) return 0 ; 
#ifdef USE_SCALAR_KMEANS_INITIALIZATION 

if ( 'blsStrongSpot && (lNuinIter>0 || acceptKM) ) break; 



#else 
#endif 



if ( IblsStrongSpot && lNumIter>0) break; 



pBinPrev=pBin ; 

//Update precomp^ratio for std. devs 

precomp_ratio=l . 0 ; 

for (c = 0 ; c<lNuinColors ; C++ ) { 

dStdDev[0] [c] = (dStdDev[0] [ c ] ?1 . 0/dStdDev [ 0 ] [c] :0.0) ; 

dStdDev[l] [ c ] = (dStdDev [ 1 ] [c ] ?1 . 0/dStdDev [ 1 ] [c] :0,0) ; 

precomp_ratio*= (dStdDev[l] [c] ?dStdDev[l] [c] :1.0) /(dStdDev[0] [c] ?dStdD 
ev[0] [c] :1.0) ; 

} 

precomp_ratio=2 . O*log (precomp^ratio) ; 
//Classify pixels 

for (p=0 , y=lYStart ; y<=lYStop; y++ ) { 

for (x=lXStart ;x<=lXStop;x++,p++) { 
metric [1 1 =pPosProb [p] ; 
metric [0]=1.0-metric [1] ; 

if (metric [1]==1) pBin[p]=l; 
else if (metric [1] ==0) pBin[p]=0; 
else { 

metric [1] =2 . O*log(metric [1] /metric [0] ) ; 
metric [0] =0; 

for (c=0 ; c<lNiamColors ;c++) { 

t= (double) ppPixelData [c] [p] ; 

metric [0] +=log_gaus si an_prob( (double) ppPixelData[c] [p] ,dAve[0] [c] ,dSt 
dDev [ 0 ] [ c ] ) ; 

metric [1] +=log_gaussian_prob ( (double) ppPixelData [c ] [p] , dAve [1] [c] , dSt 
dDev[l] [c] ) ; 

} 

pBin[p] = (metric [1] +precomp_ratio>metric [0] ) ; 

} 

} 

} 

//Convergence test 

for (bIsConverged=true,p=0;p<lNumPixels && bIsConverged;p++) 

bIsConverged= (pBinPrev[p] ==pBin[p] ) ; 
lNumIter++; 
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DoHoleFill (IXStart, IXStop, lYStart, lYStop, pBin) ; 

//Compute number of fg pixels 
^ for (ntO] =n[l] =0 ,p=0 ;p<lN\imPixels;p++) n[pBin[p] ] ++ ; 

#ifdef UPDATE_CENTROID 

//Compute centroid of fg region 
long sum_x=0, sum_y=0; 
for (p=0,y=lyStart;y<=lYStop;y++) { 
10 for (x=lXStart;x<=lXStop;x++,p++) { 

if (pBin[p3) { 

sum_x+ =x ; sum v+=v ; 

} 

} 

15 } 

if (n[l]) { 

dCentroid_x = ( (double) sum__x/n[l] ) ; 

dCentroid__y = ( (double) sum_y/n[l] ) ; 

} 

20 #endif 



//Final check for no spots 
nom_rad=sqrt ( (double) n [1 ] /3 . 14159) ; 

GetlnitialClassif icationNomRaddXStart , IXStop, lYStart, lYStop, 
25 dCentroid__x, 
dCentroid v, nom_rad, pBin) ; 

dSum[0] . Initialize ( zero) ; 
dSum[l] . Initialize { zero) ; 
dSuinSq[0] . Initialize ( zero) ; 
30 dSumSq [ 1] . Initialize ( zero) ; 

for (c=0;c<lNumColors;c++) { 

for (n[0] =n[l] =p=0;p<lNuinPixels;p++) { 
n[pBin[p] ] ++; 

dSum[pBin[p] ] [c] += (double) ppPixelData [c] [p] ; 

35 

dSumSq[pBin[p] ] [c] += (double) ppPixelData [c ] [p] * (double) ppPixelData[c] [ 

p]; 

} 

double t=dAve[0] [c] = (n [0 ] ?dSum[0 ] [c]/n[0] :0) ; 
40 t=dAve[l] [c]={n[l]?dSum[l] [c]/n[l] :0) ; 

t=dStdDev[0] [c]=(n[0]>l?sqrt (dSumSq[0] [c]- 
dAveCO] [c] *dAve[0] [c ] ) / (n [ 0 ] -1 ) : 0 ) ; 

t=dStdDev[l] [c] = (n[l]>l?sqrt (dSumSq[l] [c]- 
dAve[l] [c] *dAve[l] [c ] ) / (n [ 1 ] -1 ) : 0 ) ; 

45 } 

bool bIsSpot=f alse; 

for (c=0 ; c<lNumColors ; C++) { 

bIsSpot|=(fabs(dAve[l] [c]- 
dAve[03 [c] )>(dStdDev[l] [c] +dStdDev [0] [c] ) *ISSPOT_F ACTOR) ; 
50 } 

if ( IblsSpot) 
{ 

return 0; 

} 

55 

bIsSpotShif ted= (sqrt ( (gridCentroid_x-dCentroid_x) * (gridCentroid_x- 
dCentroid_x) + (gridCentroid v-dCentroid v) * (crridCentroid v- 
dCentroid_y) ) >GRID_MOVED_REJECT_FACTOR*nom_rad) ; 
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return n[l] ; 
} 



